Homework 5
Describe the setup and each step in your solutions with words and clearly label your final answers. Use Matlab for plotting and programming and include your code as an appendix to your problem set.
Table of Contents
Problem 1
Streeter–Phelps Model. The Streeter–Phelps equation is used in the study of water pollution as a water quality modeling tool. The model describes how dissolved oxygen (DO) decreases in a river or stream along a certain distance by degradation of biochemical oxygen demand (BOD). The equation was derived by H. W. Streeter, a sanitary engineer, and Earle B. Phelps, a consultant for the U.S. Public Health Service, in 1925, based on field data from the Ohio River. The equation is also known as the DO sag equation.
In this problem set, we consider the Streeter-Phelps model in a 2D river. The governing equations for BOD and DO includes advection and diffusion process, with additional source and sink terms.
∂ B ∂ t + u i ∂ B ∂ x i = D i ∂ 2 B ∂ x i 2 − K r B + δ ( x − x 0 , y − y 0 ) Q B , \begin{equation*}
\tag{1}
\frac{\partial B}{\partial t} + u_i \frac{\partial B}{\partial x_i} = D_i \frac{\partial^2 B}{\partial x_i^2} - K_{\text{r}} B + \delta(x - x_0, y - y_0) Q_{\text{B}},
\end{equation*} ∂ t ∂ B + u i ∂ x i ∂ B = D i ∂ x i 2 ∂ 2 B − K r B + δ ( x − x 0 , y − y 0 ) Q B , ( 1 )
∂ O ∂ t + u i ∂ O ∂ x i = D i ∂ 2 O ∂ x i 2 + K a ( O sat − O ) − K r B , \begin{equation*}
\tag{2}
\frac{\partial O}{\partial t} + u_i \frac{\partial O}{\partial x_i} = D_i \frac{\partial^2 O}{\partial x_i^2} + K_{\text{a}} (O_{\text{sat}} - O) - K_{\text{r}} B,
\end{equation*} ∂ t ∂ O + u i ∂ x i ∂ O = D i ∂ x i 2 ∂ 2 O + K a ( O sat − O ) − K r B , ( 2 )
where B B B O O O [ M L − 3 ] [\mathrm{M} \mathrm{L}^{−3}] [ M L − 3 ] K r K_{\text{r}} K r [ T − 1 ] [\mathrm{T}^{−1}] [ T − 1 ] K a K_{\text{a}} K a [ T − 1 ] [\mathrm{T}^{−1}] [ T − 1 ] K a ( O sat − O ) K_{\text{a}} (O_{\text{sat}} −O) K a ( O sat − O ) O sat O_{\text{sat}} O sat Q B Q_{\text{B}} Q B [ M T − 1 ] [\mathrm{M} \mathrm{T}^{−1}] [ M T − 1 ] δ ( x − x 0 , y − y 0 ) \delta(x −x_0, y −y_0) δ ( x − x 0 , y − y 0 ) i = 1 , 2 i = 1, 2 i = 1 , 2 x x x y y y
u i ∂ B ∂ x i ≡ u 1 ∂ B ∂ x 1 + u 2 ∂ B ∂ x 2 = u ∂ B ∂ x + v ∂ B ∂ y . \begin{equation*}
u_i \frac{\partial B}{\partial x_i} \equiv u_1 \frac{\partial B}{\partial x_1} + u_2 \frac{\partial B}{\partial x_2} = u \frac{\partial B}{\partial x} + v \frac{\partial B}{\partial y}.
\end{equation*} u i ∂ x i ∂ B ≡ u 1 ∂ x 1 ∂ B + u 2 ∂ x 2 ∂ B = u ∂ x ∂ B + v ∂ y ∂ B .
1. Discretize the above two equations using the forward in time, centered in space (FTCS) scheme. Set the proper lateral boundary conditions for the problem (HINT: use solid wall boundaries in the cross flow y y y
2. Perform von Neumann stability analysis for the discretized equations. Determine the proper proper Courant number Cr = u Δ t / Δ x , r = D Δ t / Δ x 2 \text{Cr} = u \Delta t / \Delta x, \, r = D \Delta t / {\Delta x}^2 Cr = u Δ t /Δ x , r = D Δ t / Δ x 2 Pe = u Δ x / D \text{Pe} = u \Delta x / D Pe = u Δ x / D δ ( x − x 0 , y − y 0 ) Q B \delta(x −x_0, y − y_0) Q_{\text{B}} δ ( x − x 0 , y − y 0 ) Q B K a O sat − K r B K_{\text{a}} O_{\text{sat}} − K_{\text{r}} B K a O sat − K r B D x = D y D_x = D_y D x = D y
3. Perform modified equation analysis, and discuss the truncation errors.
4. Write a Matlab script to solve the problem. Use the following assumptions:
The length and width of the river is L x = 150 m , L y = 30 m L_x = 150 \, \text{m}, \, L_y = 30 \, \text{m} L x = 150 m , L y = 30 m
The flow is assumed to be in the x x x u = 0.4 m s − 1 , v = 0 u = 0.4 \, \text{m} \, \text{s}^{−1}, v = 0 u = 0.4 m s − 1 , v = 0
The diffusion coefficients D x = D y = 0.5 m 2 s − 1 D_x = D_y = 0.5 \, \text{m}^2 \, \text{s}^{−1} D x = D y = 0.5 m 2 s − 1
The reaction rate K r = 0.01 s − 1 K_{\text{r}} = 0.01 \, \text{s}^{−1} K r = 0.01 s − 1
The mass transfer coefficient for oxygen K a = 0.02 s − 1 K_{\text{a}} = 0.02 \, \text{s}^{−1} K a = 0.02 s − 1
The saturation oxygen content O sat = 8 mg L − 1 O_{\text{sat}} = 8 \, \text{mg} \; \text{L}^{−1} O sat = 8 mg L − 1
The BOD discharge rate Q B = 70 mg s − 1 Q_{\text{B}} = 70 \, \text{mg} \; \text{s}^{−1} Q B = 70 mg s − 1
Place the point source Q B Q_{\text{B}} Q B x = 10 m , y = 15 m x = 10 \, \text{m}, \, y = 15 \, \text{m} x = 10 m , y = 15 m
Use Δ x = Δ y = 1 m \Delta x = \Delta y = 1 \, \text{m} Δ x = Δ y = 1 m
5. Plot a x , y x, y x , y ∂ / ∂ t = 0 \partial / \partial t = 0 ∂ / ∂ t = 0 contour or contourf). Plot the DO sag curve and BOD curve along the center line y = 15 m y =15 \, \text{m} y = 15 m B ( x , y = 15 m ) B(x, y = 15 \, \text{m}) B ( x , y = 15 m ) O ( x , y = 15 m ) O(x, y = 15 \, \text{m}) O ( x , y = 15 m )
6. According to Environmental Protection Agency (EPA) requirement, the minimum steady state DO in a river must be above 6 mg L − 1 6 \, \text{mg} \; \text{L}^{−1} 6 mg L − 1 70 mg s − 1 70 \, \text{mg} \; \text{s}^{−1} 70 mg s − 1
7. Suggest 3 or more ways to make your code more realisitic.
8. Suggest 3 or more ways to speed up your code.
本章讨论二维 Streeter-Phelps 模型 (1)(2) 的数值解. 用 forward in time, centered in space (FTCS) 方案, 将 Eq. (1) 离散为
1 τ ( B i , j n + 1 − B i , j n ) + u i , j n 2 h x ( B i + 1 , j n − B i − 1 , j n ) + v i , j n 2 h y ( B i , j + 1 n − B i , j − 1 n ) = D x h x 2 ( B i + 1 , j n − 2 B i , j n + B i − 1 , j n ) + D y h y 2 ( B i , j + 1 n − 2 B i , j n + B i , j − 1 n ) − K r B i , j n + Q B h x h y I i , j ( x 0 , y 0 ) , \begin{equation*}
\tag{3a}
\begin{aligned}
& \quad \; \frac{1}{\tau} \left( B^{n+1}_{i,j} - B^{n}_{i,j} \right) + \frac{u^{n}_{i, j}}{2 h_x} \left( B^{n}_{i+1,j} - B^{n}_{i-1,j} \right) + \frac{v^n_{i,j}}{2 h_y} \left( B^{n}_{i,j+1} - B^{n}_{i,j-1} \right) \\
& = \frac{D_x}{h^2_x} \left( B^{n}_{i+1,j} - 2 B^{n}_{i,j} + B^{n}_{i-1,j} \right) + \frac{D_y}{h^2_y} \left( B^{n}_{i,j+1} - 2 B^{n}_{i,j} + B^{n}_{i,j-1} \right) - K_{\text{r}} B^n_{i,j} + \frac{Q_\text{B}}{h_x h_y} \mathrm{I}_{i,j}(x_0, y_0),
\end{aligned}
\end{equation*} τ 1 ( B i , j n + 1 − B i , j n ) + 2 h x u i , j n ( B i + 1 , j n − B i − 1 , j n ) + 2 h y v i , j n ( B i , j + 1 n − B i , j − 1 n ) = h x 2 D x ( B i + 1 , j n − 2 B i , j n + B i − 1 , j n ) + h y 2 D y ( B i , j + 1 n − 2 B i , j n + B i , j − 1 n ) − K r B i , j n + h x h y Q B I i , j ( x 0 , y 0 ) , ( 3a )
where
n = 0 , 1 , ⋯ , N , i = 0 , 1 , ⋯ I , j = 0 , 1 , ⋯ , J , n = 0, 1, \cdots, N, \quad i = 0, 1, \cdots I, \quad j = 0, 1, \cdots, J, n = 0 , 1 , ⋯ , N , i = 0 , 1 , ⋯ I , j = 0 , 1 , ⋯ , J ,
示性函数
I i , j ( x 0 , y 0 ) : = { 1 , 源点 ( x 0 , y 0 ) 落在 ( i , j ) 格点所代表的那个有限体积元 , 0 , otherwise . \begin{equation*}
\tag{4}
\mathrm{I}_{i,j}(x_0, y_0) := \left\{
\begin{aligned}
& 1, \quad 源点 (x_0, y_0) \text{ 落在 $(i, j)$ 格点所代表的那个有限体积元}, \\
& 0, \quad \text{otherwise}.
\end{aligned}
\right.
\end{equation*} I i , j ( x 0 , y 0 ) := { 1 , 源点 ( x 0 , y 0 ) 落在 ( i , j ) 格点所代表的那个有限体积元 , 0 , otherwise . ( 4 )
可以设置多个源点. 将 Eq. (2) 离散为
1 τ ( O i , j n + 1 − O i , j n ) + u i , j n 2 h x ( O i + 1 , j n − O i − 1 , j n ) + v i , j n 2 h y ( O i , j + 1 n − O i , j − 1 n ) = D x h x 2 ( O i + 1 , j n − 2 O i , j n + O i − 1 , j n ) + D y h y 2 ( O i , j + 1 n − 2 O i , j n + O i , j − 1 n ) + K a ( O sat − O i , j n ) − K r B i , j n . \begin{equation*}
\tag{3b}
\begin{aligned}
& \quad \; \frac{1}{\tau} \left( O^{n+1}_{i,j} - O^{n}_{i,j} \right) + \frac{u^{n}_{i, j}}{2 h_x} \left( O^{n}_{i+1,j} - O^{n}_{i-1,j} \right) + \frac{v^n_{i,j}}{2 h_y} \left( O^{n}_{i,j+1} - O^{n}_{i,j-1} \right) \\
& = \frac{D_x}{h^2_x} \left( O^{n}_{i+1,j} - 2 O^{n}_{i,j} + O^{n}_{i-1,j} \right) + \frac{D_y}{h^2_y} \left( O^{n}_{i,j+1} - 2 O^{n}_{i,j} + O^{n}_{i,j-1} \right) + K_{\text{a}} \left( O_{\text{sat}} - O^n_{i,j} \right) - K_{\text{r}} B^n_{i,j}.
\end{aligned}
\end{equation*} τ 1 ( O i , j n + 1 − O i , j n ) + 2 h x u i , j n ( O i + 1 , j n − O i − 1 , j n ) + 2 h y v i , j n ( O i , j + 1 n − O i , j − 1 n ) = h x 2 D x ( O i + 1 , j n − 2 O i , j n + O i − 1 , j n ) + h y 2 D y ( O i , j + 1 n − 2 O i , j n + O i , j − 1 n ) + K a ( O sat − O i , j n ) − K r B i , j n . ( 3b )
格式 (3) 中的速度场 u n , v n u^n, v^n u n , v n
Problem 1.1
我们要在二维有界区域
t ≥ 0 , ( x , y ) ∈ Ω : = [ 0 , L x ] × [ 0 , L y ] \begin{equation*}
\tag{5}
t \ge 0, \quad (x, y) \in \Omega := [0, L_x] \times [0, L_y]
\end{equation*} t ≥ 0 , ( x , y ) ∈ Ω := [ 0 , L x ] × [ 0 , L y ] ( 5 )
上求解模型 (1)(2), 为此, 须为模型设置合适的初始条件和边界条件, 这里取为
{ ( I.C. ) F ∣ t = 0 = F 0 ( x , y ) , ( x , y ) ∈ Ω , ( B.C. ) { ( Dirichlet ) F ∣ x = 0 = F 0 ( y , t ) , ( Neumann ) F x ∣ x = L x = F 1 ( y , t ) , ( solid wall ) ( n ⋅ u ) y = 0 , L y = 0 , ( solid wall 2 ) [ n ⋅ ( − κ ∇ F ) ] y = 0 , L y = 0. \begin{equation*}
\tag{6}
\left\{
\begin{aligned}
& (\text{I.C.}) \quad
\begin{aligned}
% \tag{6}
\left. F \right|_{t = 0} = F^0(x, y), \quad (x,y) \in \Omega,
\end{aligned} \\
& (\text{B.C.}) \quad
\left\{
\begin{aligned}
% \tag{7}
& (\text{Dirichlet}) & \quad \left. F \right|_{x = 0} = F_0(y, t), \\
& (\text{Neumann}) & \quad \left. F_x \right|_{x = L_x} = F_1(y, t), \\
& (\text{solid wall}) & \quad \left( \boldsymbol{n} \cdot \boldsymbol{u} \right)_{y = 0, \, L_y} = 0, \\
& ({\color{blue}{\text{solid wall 2}}}) & \quad \left[ \boldsymbol{n} \cdot (- \kappa \nabla F) \right]_{y = 0, \, L_y} = 0.
\end{aligned}
\right. \\
\end{aligned}
\right.
\end{equation*} ⎩ ⎨ ⎧ ( I.C. ) F ∣ t = 0 = F 0 ( x , y ) , ( x , y ) ∈ Ω , ( B.C. ) ⎩ ⎨ ⎧ ( Dirichlet ) ( Neumann ) ( solid wall ) ( solid wall 2 ) F ∣ x = 0 = F 0 ( y , t ) , F x ∣ x = L x = F 1 ( y , t ) , ( n ⋅ u ) y = 0 , L y = 0 , [ n ⋅ ( − κ ∇ F ) ] y = 0 , L y = 0. ( 6 )
Eq. (6) 中, F F F B B B O O O
相应地, 为格式 (3) 引入由
F − 1 , j n + F 1 , j n 2 = F 0 ( t n , y j ) , F 1 , j n − F − 1 , j n 2 h x = F 1 ( t n , y j ) , F i , − 1 n = 3 F i , 0 n − 3 F i , 1 n + F i , 2 n , F i , J + 1 n = 3 F i , J n − 3 F i , J − 1 n + F i , J − 2 n , \begin{align*}
\tag{7a}
& \frac{F^n_{-1, j} + F^n_{1, j}}{2} = F_0(t_n, y_j), \\
\tag{7b}
& \frac{F^n_{1, j} - F^n_{-1, j}}{2 h_x} = F_1(t_n, y_j), \\
\tag{8a}
& \color{gray}{F^n_{i, -1} = 3 F^n_{i, 0} - 3F^n_{i, 1} + F^n_{i, 2}}, \\
\tag{8b}
& \color{gray}{F^n_{i, J+1} = 3 F^n_{i, J} - 3 F^n_{i, J-1} + F^n_{i, J-2}},
\end{align*} 2 F − 1 , j n + F 1 , j n = F 0 ( t n , y j ) , 2 h x F 1 , j n − F − 1 , j n = F 1 ( t n , y j ) , F i , − 1 n = 3 F i , 0 n − 3 F i , 1 n + F i , 2 n , F i , J + 1 n = 3 F i , J n − 3 F i , J − 1 n + F i , J − 2 n , ( 7a ) ( 7b ) ( 8a ) ( 8b )
决定的 halo points F − 1 , j n , F I + 1 , j n , F i , − 1 n , F i , J + 1 n F^n_{-1, j}, F^n_{I + 1, j}, F^n_{i, -1}, F^n_{i, J + 1} F − 1 , j n , F I + 1 , j n , F i , − 1 n , F i , J + 1 n 其中, Eq. (8) 的设计考虑是使格式 (3) 在网格边界 j = 0 , J j = 0, J j = 0 , J O ( h y ) \mathcal{O}(h_y) O ( h y ) y = 0 , L y y = 0, L_y y = 0 , L y F i , − 1 n , F i , J + 1 n F^n_{i, -1}, F^n_{i, J+1} F i , − 1 n , F i , J + 1 n v i , 0 n = v i , J n = 0 v^n_{i, 0} = v^n_{i, J} = 0 v i , 0 n = v i , J n = 0 F i , − 1 n , F i , J + 1 n F^n_{i, -1}, F^n_{i, J+1} F i , − 1 n , F i , J + 1 n y = 0 , L y y = 0, L_y y = 0 , L y v v v exactly zero) 的离散 . 一个重要问题是, 在物理上, solid wall 边界对 Fick 扩散的约束是怎样的? 若是梯度扩散通量为零, 即 n ⋅ ( − κ ∇ F ) = 0 \color{black}{ \boldsymbol{n} \cdot (- \kappa \nabla F) = 0} n ⋅ ( − κ ∇ F ) = 0 ∂ F / ∂ n = 0 \color{black}{\partial F / \partial n = 0} ∂ F / ∂ n = 0 发现 ChatGPT 支持 此观点, 且有符合常识的物理解释, 所以我将 (8) 改为
F i , − 1 n = F i , 1 n , F i , J + 1 n = F i , J − 1 n . \begin{equation*}
\tag{$8'$}
F^n_{i, -1} = F^n_{i, 1}, \quad F^n_{i, J+1} = F^n_{i, J-1}.
\end{equation*} F i , − 1 n = F i , 1 n , F i , J + 1 n = F i , J − 1 n . ( 8 ′ )
但, 若我坚持用 (8), 则会引入多大的误差呢? 毕竟, Eq. (8), 尽管不蕴含 (8'), 似乎也是在合理前提下 (∂ 2 F / ∂ n 2 \partial^2 F / \partial n^2 ∂ 2 F / ∂ n 2 导出的?
引入 halo points (7)(8') 后, 在速度场 u i , j n , v i , j n u^n_{i,j}, v^n_{i,j} u i , j n , v i , j n
B i , j 0 = B 0 ( x i , y j ) , O i , j 0 = O 0 ( x i , y j ) , x i = i h x , y j = j h y , \begin{equation*}
\tag{9}
B^0_{i,j} = B^0(x_i, y_j), \quad O^0_{i,j} = O^0(x_i, y_j), \quad x_i = i h_x, \, y_j = j h_y,
\end{equation*} B i , j 0 = B 0 ( x i , y j ) , O i , j 0 = O 0 ( x i , y j ) , x i = i h x , y j = j h y , ( 9 )
就可驱动格式 (3) 向前积分. 以后, 简单起见, 将速度场取为常向量, 即认为 u n ≡ u , v n ≡ v u^n \equiv u, \, v^n \equiv v u n ≡ u , v n ≡ v u , v u, v u , v
Problem 1.2
本节对格式 (3) 作 von Neumann 稳定性分析. 为简化问题, I ignored the BOD discharge term δ ( x − x 0 , y − y 0 ) Q B \delta(x −x_0, y − y_0) Q_{\text{B}} δ ( x − x 0 , y − y 0 ) Q B K a O sat − K r B K_{\text{a}} O_{\text{sat}} − K_{\text{r}} B K a O sat − K r B
G = ( 1 − 2 r x − 2 r y − K τ ) + 2 r x cos ( ξ x h x ) − i C r x sin ( ξ x h x ) + 2 r y cos ( ξ y h y ) − i C r y sin ( ξ y h y ) , \begin{equation*}
\tag{10}
\begin{aligned}
G & = (1 - 2 r_x - 2 r_y - K \tau) \\
& + 2 r_x \cos{(\xi_x h_x)} - \mathrm{i} \, \mathrm{Cr}_x \, \sin{(\xi_x h_x)} \\
& + 2 r_y \cos{(\xi_y h_y)} - \mathrm{i} \, \mathrm{Cr}_y \, \sin{(\xi_y h_y)},
\end{aligned}
\end{equation*} G = ( 1 − 2 r x − 2 r y − K τ ) + 2 r x cos ( ξ x h x ) − i Cr x sin ( ξ x h x ) + 2 r y cos ( ξ y h y ) − i Cr y sin ( ξ y h y ) , ( 10 )
其中 K K K K r K_{\mathrm{r}} K r K a K_{\mathrm{a}} K a
C r i : = u i τ / h i 2 , r i : = D τ / h i 2 , P e , i : = u i h / D , \begin{equation*}
\mathrm{Cr}_i := u_i \tau / h_i^2, \; r_i := D \tau / h_i^2, \; P_{\mathrm{e}, i} := u_i h / D,
\end{equation*} Cr i := u i τ / h i 2 , r i := D τ / h i 2 , P e , i := u i h / D ,
下标 i = 1 , 2 , 3 i = 1, 2, 3 i = 1 , 2 , 3 x , y , z x, y, z x , y , z τ , h x , h y \tau, h_x, h_y τ , h x , h y Hw 3 的 Eq. (1.22))
∣ G ∣ ≤ 1 , ∀ ∣ ξ x h x ∣ ≤ π , ∣ ξ y h y ∣ ≤ π \begin{equation*}
\tag{11}
|G| \le 1, \quad \forall |\xi_x h_x| \le \pi, \, |\xi_y h_y| \le \pi
\end{equation*} ∣ G ∣ ≤ 1 , ∀∣ ξ x h x ∣ ≤ π , ∣ ξ y h y ∣ ≤ π ( 11 )
被满足.
设 A A A A 1 , A 2 ⊂ A A_1, A_2 \subset A A 1 , A 2 ⊂ A
A 1 + A 2 : = { a 1 + a 2 ∣ a i ∈ A i } . \begin{equation*}
\tag{12}
A_1 + A_2 := \left\{ \left. a_1 + a_2 \; \right| \; a_i \in A_i \right\}.
\end{equation*} A 1 + A 2 := { a 1 + a 2 ∣ a i ∈ A i } . ( 12 )
记平面直角坐标系上的椭圆
C ( a , b ) : = { ( x , y ) ∣ x 2 a 2 + y 2 b 2 = 1 } , a , b > 0. \begin{equation*}
\tag{13}
C(a, b) := \left\{(x, y) \; \left| \; \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 \right. \right\}, \quad a, b > 0.
\end{equation*} C ( a , b ) := { ( x , y ) a 2 x 2 + b 2 y 2 = 1 } , a , b > 0. ( 13 )
当 ( ξ x h x , ξ y h y ) (\xi_x h_x, \xi_y h_y) ( ξ x h x , ξ y h y ) [ − π , π ) 2 [-\pi, \pi)^2 [ − π , π ) 2 G = x + i y G = x + \mathrm{i} y G = x + i y
z 0 + C ( 2 r x , ∣ C r x ∣ ) + C ( 2 r y , ∣ C r y ∣ ) , \begin{equation*}
\tag{14}
z_0 + C(2 r_x, |\mathrm{Cr}_x|) + C(2 r_y, |\mathrm{Cr}_y|),
\end{equation*} z 0 + C ( 2 r x , ∣ Cr x ∣ ) + C ( 2 r y , ∣ Cr y ∣ ) , ( 14 )
其中
z 0 : = 1 − 2 r x − 2 r y − K τ \begin{equation*}
\tag{15}
z_0 := 1 - 2 r_x - 2 r_y - K \tau
\end{equation*} z 0 := 1 − 2 r x − 2 r y − K τ ( 15 )
在实轴上. 由 (14) 看出, 命题 (11) 的一个充分条件是
∣ z 0 ∣ + max { 2 r x , ∣ C r x ∣ } + max { 2 r y , ∣ C r y ∣ } ≤ 1 , \begin{equation*}
\tag{16}
|z_0| + \max{\{2 r_x, |\mathrm{Cr}_x|\}} + \max{\{2 r_y, |\mathrm{Cr}_y|\}} \le 1,
\end{equation*} ∣ z 0 ∣ + max { 2 r x , ∣ Cr x ∣ } + max { 2 r y , ∣ Cr y ∣ } ≤ 1 , ( 16 )
一个必要条件是
∣ z 0 ∣ + 2 r x + 2 r y ≤ 1. \begin{equation*}
\tag{17}
|z_0| + 2 r_x + 2 r_y \le 1.
\end{equation*} ∣ z 0 ∣ + 2 r x + 2 r y ≤ 1. ( 17 )
故在
2 r x ≥ ∣ C r x ∣ , 2 r y ≥ ∣ C r y ∣ \begin{equation*}
\tag{18}
2 r_x \ge |\mathrm{Cr}_x|, \, 2 r_y \ge |\mathrm{Cr}_y|
\end{equation*} 2 r x ≥ ∣ Cr x ∣ , 2 r y ≥ ∣ Cr y ∣ ( 18 )
的大前提下, 命题 (11) 的充分条件 (16) 即是必要条件 (17), 从而是充要的. 对于 (18) 以外的情形, (16) 是充分但不必要的. 简单起见, 我们用命题 (16) 的、从而也是命题 (11) 的一个充分不必要 条件,
∣ Cr y ∣ ≤ ∣ Cr x ∣ ≤ 2 r x = 2 r y ≤ ( 1 − K τ ) / 2 , \begin{equation*}
\tag{19}
|\text{Cr}_y| \le |\text{Cr}_x| \le 2 r_x = 2 r_y \le (1 - K \tau) / 2,
\end{equation*} ∣ Cr y ∣ ≤ ∣ Cr x ∣ ≤ 2 r x = 2 r y ≤ ( 1 − K τ ) /2 , ( 19 )
来决定网格参数 h x , h y , τ h_x, h_y, \tau h x , h y , τ v = 0 , D x = D y = D v = 0, \, D_x = D_y = D v = 0 , D x = D y = D
h x = h y = : h ≤ 2 D / ∣ u ∣ , τ ≤ ( 4 D / h 2 + max { K a , K r } ) − 1 . \begin{equation*}
h_x = h_y =: h \le {2 D / |u|}, \quad \tau \le \left(4 D / h^2 + \max{\{ K_{\mathrm{a}}, K_{\mathrm{r}} \}} \right)^{-1}.
\end{equation*} h x = h y =: h ≤ 2 D /∣ u ∣ , τ ≤ ( 4 D / h 2 + max { K a , K r } ) − 1 .
一般地, 要写出 (11) 的充要条件, 或许可先用 Lagrange multipliers (2019 , 2022 ) 方法写出 max ∣ ξ h ∣ ≤ π ∣ G ∣ \max_{|\xi h| \le \pi}|G| max ∣ ξ h ∣ ≤ π ∣ G ∣ G G G MMA ) 软件算得 (11) 的一个充要条件是
MMA 似乎算不出来? 是我给的命令不好吗? \begin{equation*}
\tag{20}
\text{MMA 似乎算不出来? 是我给的命令不好吗?}
\end{equation*} MMA 似乎算不出来 ? 是我给的命令不好吗 ? ( 20 )
G = 1 - 2 rx - 2 ry - k tau + 2 rx cx - I crx sx + 2 ry cy - I crx sy
ForAll[ { cx, sx, cy, sy} , cx^ 2 + sx^ 2 == 1 & & cy^ 2 + sy^ 2 == 1 , Abs [ G] <= 1 ]
若作简化
Cr y = 0 , r x = r y = r , K τ = 0 , \begin{equation*}
\tag{21}
\text{Cr}_y = 0, \; r_x = r_y = r, \; K\tau = 0,
\end{equation*} Cr y = 0 , r x = r y = r , K τ = 0 , ( 21 )
也就是将命题 (21) 添加到大前提, 则用 MMA 显式化命题 (11)
Resolve[ ForAll[ { cx, sx, cy} , cx^ 2 + sx^ 2 == 1 & & - 1 <= cy <= 1 , ( crx* sx) ^ 2 + ( 2 * cx* r + 2 * cy* r - 2 * r - 2 * r + 1 ) ^ 2 <= 1 ] , Reals]
FullSimplify[ %, r >= 0 ]
将得到形式很复杂的复合命题
crx 2 + 4 ( r − 1 ) r ≤ 0 ∧ 4 r ≤ 1 ∧ ( crx = 2 r ∨ r = 0 ∨ crx + 2 r = 0 ∨ ( ( ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≤ 0 ∨ ( r 2 ( 2 r − 1 ) ( 4 r 2 − crx 2 ) ≤ 0 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ) ∨ 6 r < 1 ∨ 2 r 2 ≥ r ∨ r 2 ≤ 0 ) ∧ ( crx ≠ 0 ∨ crx 2 r ( 4 r − 1 ) ≥ 0 ∨ crx 2 r ( 4 r − 1 ) ( crx 4 − 4 crx 2 r ( r + 1 ) + 4 r 2 ) ≥ 0 ∨ ( r 2 ( crx 10 − 4 crx 8 ( r ( 9 r − 4 ) + 1 ) + 8 crx 6 r 2 ( 62 r 2 − 48 r + 11 ) + 16 crx 4 r 2 ( 4 r ( r ( 3 r ( 3 r − 4 ) + 10 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≥ 0 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ) ∨ ( r 2 ( crx 10 + 4 crx 8 r ( 7 r − 4 ) − 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≤ 0 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 4 + 4 crx 2 r ( 3 r − 2 ) + 8 r 2 ( 8 r 2 − 4 r + 1 ) ≠ 0 ) ∧ ( crx 10 + 4 crx 8 r ( 7 r − 4 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) + 64 r 6 ≠ 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) ∨ crx 2 r 2 + 4 r 4 ≤ 0 ∨ ( crx 2 + 4 r 2 ≠ 0 ∧ crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ≠ 0 ) ) ) ∨ ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 ≠ 0 ∧ crx 2 r ( 4 r − 1 ) ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 ) ≤ 0 ) ∨ ( crx 2 r ( 4 r − 1 ) ( crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ) ≤ 0 ∧ ( r 2 ≤ 0 ∨ crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ≠ 0 ) ) ) ∧ ( crx = 0 ∨ ( crx 2 = 2 r ∧ r > 0 ) ∨ ( ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ∧ crx 8 r 2 + 16 crx 2 r 6 ≤ 4 crx 4 ( crx 2 + 1 ) r 4 ) ∨ r = 0 ∨ ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≠ 0 ∧ crx 2 r 2 ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ) ≥ 0 ∧ crx 2 r ( 6 r − 1 ) ≤ 0 ) ∨ crx 2 r ( 2 r − 1 ) ≥ 0 ∨ crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ) ∧ ( 2 r 2 ≥ r ∨ ( crx 2 − 2 r ) 2 ≤ 0 ∨ r ( 2 r − 1 ) ( crx 2 − 2 r ) > 0 ) ∧ ( ( r > 0 ∧ ( r ( 2 r − 1 ) ( crx 2 − 2 r ) < 0 ∨ ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) > 0 ∧ 6 r < 1 ) ) ) ∨ crx 2 + 4 r 2 = 0 ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) ≤ 0 ∨ crx 8 r 2 + 64 r 8 ≤ 4 crx 4 ( 4 r 6 + r 4 ) ) ) ) ∧ ( crx = 2 r ∨ r = 0 ∨ crx + 2 r = 0 ∨ ( ( crx ≠ 0 ∨ crx 2 r ( 4 r − 1 ) ≥ 0 ∨ crx 2 r ( 4 r − 1 ) ( crx 4 − 4 crx 2 r ( r + 1 ) + 4 r 2 ) ≥ 0 ∨ ( r 2 ( crx 10 − 4 crx 8 ( r ( 9 r − 4 ) + 1 ) + 8 crx 6 r 2 ( 62 r 2 − 48 r + 11 ) + 16 crx 4 r 2 ( 4 r ( r ( 3 r ( 3 r − 4 ) + 10 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≥ 0 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ) ∨ ( r 2 ( crx 10 + 4 crx 8 r ( 7 r − 4 ) − 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≤ 0 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 4 + 4 crx 2 r ( 3 r − 2 ) + 8 r 2 ( 8 r 2 − 4 r + 1 ) ≠ 0 ) ∧ ( crx 10 + 4 crx 8 r ( 7 r − 4 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) + 64 r 6 ≠ 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) ∨ crx 2 r 2 + 4 r 4 ≤ 0 ∨ ( crx 2 + 4 r 2 ≠ 0 ∧ crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ≠ 0 ) ) ) ∨ ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 ≠ 0 ∧ crx 2 r ( 4 r − 1 ) ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 ) ≤ 0 ) ∨ crx 2 r ( 4 r − 1 ) ( crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ) ≥ 0 ) ∧ ( crx = 0 ∨ ( crx 2 = 2 r ∧ r > 0 ) ∨ ( ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ∧ crx 8 r 2 + 16 crx 2 r 6 ≤ 4 crx 4 ( crx 2 + 1 ) r 4 ) ∨ ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≠ 0 ∧ crx 2 r 2 ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ) ≥ 0 ∧ crx 2 r ( 6 r − 1 ) ≤ 0 ) ∨ crx 2 r ( 2 r − 1 ) ≥ 0 ∨ crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ∨ crx 2 r 2 > 0 ) ∧ ( 6 r 2 ≤ r ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≤ 0 ∨ r ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) ) > 0 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) < 0 ) ∧ ( ( r > 0 ∧ ( r ( 2 r − 1 ) ( crx 2 − 2 r ) < 0 ∨ ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) > 0 ∧ 6 r < 1 ) ) ) ∨ crx 2 + 4 r 2 = 0 ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) ≤ 0 ∨ crx 2 r 2 + 4 r 4 > 0 ∨ crx 8 r 2 + 64 r 8 ≤ 4 crx 4 ( 4 r 6 + r 4 ) ) ) ) ∧ ( r = 0 ∨ crx 2 ≥ 4 r 2 ∨ ( ( crx = 0 ∨ crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ∨ ( ( crx 2 r ( 6 r − 1 ) ≤ 0 ∨ ( r > 0 ∧ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≤ 0 ) ) ∧ ( crx 2 r ( 6 r − 1 ) < 0 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≠ 0 ) ) ∨ crx 2 r ( 2 r − 1 ) ≥ 0 ∨ ( crx 8 r 2 + 16 crx 2 r 6 ≤ 4 crx 4 ( crx 2 + 1 ) r 4 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ∧ ( crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ( crx 2 + 4 r 2 ) 2 ≥ 0 ∨ ( crx ≠ 2 r ∧ crx 2 ≠ 2 r ∧ crx + 2 r ≠ 0 ∧ crx 2 + 2 r ≠ 0 ) ) ) ) ∧ ( 6 r = 1 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≥ 0 ∨ r ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) ) > 0 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) < 0 ) ∧ ( crx = 2 r ∨ ( 0 < r < 1 6 ∧ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) > 0 ) ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ ( crx 2 − 1 ) r 2 ( crx 4 − 16 r 4 ) ≥ 0 ) ∧ ( crx = 2 r ∨ 0 < r < 1 6 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) < 0 ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ ( crx 2 − 1 ) r 2 ( crx 4 − 16 r 4 ) ≥ 0 ) ∧ ( crx = 2 r ∨ r ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) ) < 0 ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ 0 < r < 1 6 ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) < 0 ∨ 4 crx 4 ( 4 r 6 + r 4 ) ≤ crx 8 r 2 + 64 r 8 ) ∧ ( crx = 2 r ∨ r ( 2 r − 1 ) ( crx 2 − 2 r ) < 0 ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ ( crx 2 ( 3 − 10 r ) < 2 r ( 2 r ( 8 r − 3 ) + 1 ) ∧ 6 r < 1 ) ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) < 0 ∨ 4 crx 4 ( 4 r 6 + r 4 ) ≤ crx 8 r 2 + 64 r 8 ) ∧ ( 6 r < 1 ∨ 2 r 2 ≥ r ∨ ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ∨ ( ( crx 2 + 4 r 2 ≠ 0 ∨ crx 4 + 4 crx 2 r ( 3 r − 2 ) + 8 r 2 ( 8 r 2 − 4 r + 1 ) ≠ 0 ) ∧ ( ( crx ≠ 2 r ∧ crx + 2 r ≠ 0 ) ∨ crx 2 r 2 + 4 r 4 ≤ 0 ) ) ∨ ( r 2 ( 2 r − 1 ) ( 4 r 2 − crx 2 ) ≤ 0 ∧ ( crx 2 + 4 r 2 ≠ 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) ≠ 0 ) ) ) ) ) \begin{equation*}
\tag{22}
\text{crx}^2+4 (r-1) r\leq 0\land 4 r\leq 1\land \left(\text{crx}=2 r\lor r=0\lor \text{crx}+2 r=0\lor \left(\left(\left(\text{crx}^2-1\right) r^2 \left(\text{crx}^2-4 r^2\right)\leq 0\lor \left(r^2 (2 r-1) \left(4 r^2-\text{crx}^2\right)\leq 0\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\right)\lor 6 r<1\lor 2 r^2\geq r\lor r^2\leq 0\right)\land \left(\text{crx}\neq 0\lor \text{crx}^2 r (4 r-1)\geq 0\lor \text{crx}^2 r (4 r-1) \left(\text{crx}^4-4 \text{crx}^2 r (r+1)+4 r^2\right)\geq 0\lor \left(r^2 \left(\text{crx}^{10}-4 \text{crx}^8 (r (9 r-4)+1)+8 \text{crx}^6 r^2 \left(62 r^2-48 r+11\right)+16 \text{crx}^4 r^2 (4 r (r (3 r (3 r-4)+10)-3)+1)-16 \text{crx}^2 r^4 (8 r (5 r-2)+3)+64 r^6\right)\geq 0\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\right)\lor \left(r^2 \left(\text{crx}^{10}+4 \text{crx}^8 r (7 r-4)-8 \text{crx}^6 r^2 (2 r (r+8)-7)+16 \text{crx}^4 r^2 (4 r (r (5 r (5 r-4)+11)-3)+1)-16 \text{crx}^2 r^4 (8 r (5 r-2)+3)+64 r^6\right)\leq 0\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^4+4 \text{crx}^2 r (3 r-2)+8 r^2 \left(8 r^2-4 r+1\right)\neq 0\right)\land \left(\text{crx}^{10}+4 \text{crx}^8 r (7 r-4)+16 \text{crx}^4 r^2 (4 r (r (5 r (5 r-4)+11)-3)+1)+64 r^6\neq 8 \text{crx}^6 r^2 (2 r (r+8)-7)+16 \text{crx}^2 r^4 (8 r (5 r-2)+3)\lor \text{crx}^2 r^2+4 r^4\leq 0\lor \left(\text{crx}^2+4 r^2\neq 0\land \text{crx}^4+\text{crx}^2 (4 (2-5 r) r-2)+4 r^2\neq 0\right)\right)\right)\lor \left(\text{crx}^4+4 \text{crx}^2 r (7 r-3)+4 r^2\neq 0\land \text{crx}^2 r (4 r-1) \left(\text{crx}^4+4 \text{crx}^2 r (7 r-3)+4 r^2\right)\leq 0\right)\lor \left(\text{crx}^2 r (4 r-1) \left(\text{crx}^4+\text{crx}^2 (4 (2-5 r) r-2)+4 r^2\right)\leq 0\land \left(r^2\leq 0\lor \text{crx}^4+\text{crx}^2 (4 (2-5 r) r-2)+4 r^2\neq 0\right)\right)\right)\land \left(\text{crx}=0\lor \left(\text{crx}^2=2 r\land r>0\right)\lor \left(\left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\land \text{crx}^8 r^2+16 \text{crx}^2 r^6\leq 4 \text{crx}^4 \left(\text{crx}^2+1\right) r^4\right)\lor r=0\lor \left(\text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\neq 0\land \text{crx}^2 r^2 \left(\text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\right)\geq 0\land \text{crx}^2 r (6 r-1)\leq 0\right)\lor \text{crx}^2 r (2 r-1)\geq 0\lor \text{crx}^2 \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^2-4 r^2\right)\geq 0\right)\land \left(2 r^2\geq r\lor \left(\text{crx}^2-2 r\right)^2\leq 0\lor r (2 r-1) \left(\text{crx}^2-2 r\right)>0\right)\land \left(\left(r>0\land \left(r (2 r-1) \left(\text{crx}^2-2 r\right)<0\lor \left(\text{crx}^2 (10 r-3)+2 r (2 r (8 r-3)+1)>0\land 6 r<1\right)\right)\right)\lor \text{crx}^2+4 r^2=0\lor r^2 (2 r-1) \left(16 r^4-\text{crx}^4\right)\leq 0\lor \text{crx}^8 r^2+64 r^8\leq 4 \text{crx}^4 \left(4 r^6+r^4\right)\right)\right)\right)\land \left(\text{crx}=2 r\lor r=0\lor \text{crx}+2 r=0\lor \left(\left(\text{crx}\neq 0\lor \text{crx}^2 r (4 r-1)\geq 0\lor \text{crx}^2 r (4 r-1) \left(\text{crx}^4-4 \text{crx}^2 r (r+1)+4 r^2\right)\geq 0\lor \left(r^2 \left(\text{crx}^{10}-4 \text{crx}^8 (r (9 r-4)+1)+8 \text{crx}^6 r^2 \left(62 r^2-48 r+11\right)+16 \text{crx}^4 r^2 (4 r (r (3 r (3 r-4)+10)-3)+1)-16 \text{crx}^2 r^4 (8 r (5 r-2)+3)+64 r^6\right)\geq 0\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\right)\lor \left(r^2 \left(\text{crx}^{10}+4 \text{crx}^8 r (7 r-4)-8 \text{crx}^6 r^2 (2 r (r+8)-7)+16 \text{crx}^4 r^2 (4 r (r (5 r (5 r-4)+11)-3)+1)-16 \text{crx}^2 r^4 (8 r (5 r-2)+3)+64 r^6\right)\leq 0\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^4+4 \text{crx}^2 r (3 r-2)+8 r^2 \left(8 r^2-4 r+1\right)\neq 0\right)\land \left(\text{crx}^{10}+4 \text{crx}^8 r (7 r-4)+16 \text{crx}^4 r^2 (4 r (r (5 r (5 r-4)+11)-3)+1)+64 r^6\neq 8 \text{crx}^6 r^2 (2 r (r+8)-7)+16 \text{crx}^2 r^4 (8 r (5 r-2)+3)\lor \text{crx}^2 r^2+4 r^4\leq 0\lor \left(\text{crx}^2+4 r^2\neq 0\land \text{crx}^4+\text{crx}^2 (4 (2-5 r) r-2)+4 r^2\neq 0\right)\right)\right)\lor \left(\text{crx}^4+4 \text{crx}^2 r (7 r-3)+4 r^2\neq 0\land \text{crx}^2 r (4 r-1) \left(\text{crx}^4+4 \text{crx}^2 r (7 r-3)+4 r^2\right)\leq 0\right)\lor \text{crx}^2 r (4 r-1) \left(\text{crx}^4+\text{crx}^2 (4 (2-5 r) r-2)+4 r^2\right)\geq 0\right)\land \left(\text{crx}=0\lor \left(\text{crx}^2=2 r\land r>0\right)\lor \left(\left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\land \text{crx}^8 r^2+16 \text{crx}^2 r^6\leq 4 \text{crx}^4 \left(\text{crx}^2+1\right) r^4\right)\lor \left(\text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\neq 0\land \text{crx}^2 r^2 \left(\text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\right)\geq 0\land \text{crx}^2 r (6 r-1)\leq 0\right)\lor \text{crx}^2 r (2 r-1)\geq 0\lor \text{crx}^2 \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^2-4 r^2\right)\geq 0\lor \text{crx}^2 r^2>0\right)\land \left(6 r^2\leq r\lor \text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\leq 0\lor r \left(\text{crx}^2 (10 r-3)+2 r (2 r (8 r-3)+1)\right)>0\lor \text{crx}^4+4 \text{crx}^2 r (8 r-3)+8 r^2 (6 r (3 r-1)+1)<0\right)\land \left(\left(r>0\land \left(r (2 r-1) \left(\text{crx}^2-2 r\right)<0\lor \left(\text{crx}^2 (10 r-3)+2 r (2 r (8 r-3)+1)>0\land 6 r<1\right)\right)\right)\lor \text{crx}^2+4 r^2=0\lor r^2 (2 r-1) \left(16 r^4-\text{crx}^4\right)\leq 0\lor \text{crx}^2 r^2+4 r^4>0\lor \text{crx}^8 r^2+64 r^8\leq 4 \text{crx}^4 \left(4 r^6+r^4\right)\right)\right)\right)\land \left(r=0\lor \text{crx}^2\geq 4 r^2\lor \left(\left(\text{crx}=0\lor \text{crx}^2 \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^2-4 r^2\right)\geq 0\lor \left(\left(\text{crx}^2 r (6 r-1)\leq 0\lor \left(r>0\land \text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\leq 0\right)\right)\land \left(\text{crx}^2 r (6 r-1)<0\lor \text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\neq 0\right)\right)\lor \text{crx}^2 r (2 r-1)\geq 0\lor \left(\text{crx}^8 r^2+16 \text{crx}^2 r^6\leq 4 \text{crx}^4 \left(\text{crx}^2+1\right) r^4\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\land \left(\text{crx}^2 \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^2-4 r^2\right) \left(\text{crx}^2+4 r^2\right)^2\geq 0\lor \left(\text{crx}\neq 2 r\land \text{crx}^2\neq 2 r\land \text{crx}+2 r\neq 0\land \text{crx}^2+2 r\neq 0\right)\right)\right)\right)\land \left(6 r=1\lor \text{crx}^4+4 \text{crx}^2 r (8 r-3)+4 r^2\geq 0\lor r \left(\text{crx}^2 (10 r-3)+2 r (2 r (8 r-3)+1)\right)>0\lor \text{crx}^4+4 \text{crx}^2 r (8 r-3)+8 r^2 (6 r (3 r-1)+1)<0\right)\land \left(\text{crx}=2 r\lor \left(0<r<\frac{1}{6}\land \text{crx}^4+4 \text{crx}^2 r (8 r-3)+8 r^2 (6 r (3 r-1)+1)>0\right)\lor \text{crx}+2 r=0\lor \text{crx}^2+4 r^2=0\lor \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^4-16 r^4\right)\geq 0\right)\land \left(\text{crx}=2 r\lor 0<r<\frac{1}{6}\lor \text{crx}^4+4 \text{crx}^2 r (8 r-3)+8 r^2 (6 r (3 r-1)+1)<0\lor \text{crx}+2 r=0\lor \text{crx}^2+4 r^2=0\lor \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^4-16 r^4\right)\geq 0\right)\land \left(\text{crx}=2 r\lor r \left(\text{crx}^2 (10 r-3)+2 r (2 r (8 r-3)+1)\right)<0\lor \text{crx}+2 r=0\lor \text{crx}^2+4 r^2=0\lor 0<r<\frac{1}{6}\lor r^2 (2 r-1) \left(16 r^4-\text{crx}^4\right)<0\lor 4 \text{crx}^4 \left(4 r^6+r^4\right)\leq \text{crx}^8 r^2+64 r^8\right)\land \left(\text{crx}=2 r\lor r (2 r-1) \left(\text{crx}^2-2 r\right)<0\lor \text{crx}+2 r=0\lor \text{crx}^2+4 r^2=0\lor \left(\text{crx}^2 (3-10 r)<2 r (2 r (8 r-3)+1)\land 6 r<1\right)\lor r^2 (2 r-1) \left(16 r^4-\text{crx}^4\right)<0\lor 4 \text{crx}^4 \left(4 r^6+r^4\right)\leq \text{crx}^8 r^2+64 r^8\right)\land \left(6 r<1\lor 2 r^2\geq r\lor \left(\text{crx}^2-1\right) r^2 \left(\text{crx}^2-4 r^2\right)\geq 0\lor \left(\left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^4+4 \text{crx}^2 r (3 r-2)+8 r^2 \left(8 r^2-4 r+1\right)\neq 0\right)\land \left((\text{crx}\neq 2 r\land \text{crx}+2 r\neq 0)\lor \text{crx}^2 r^2+4 r^4\leq 0\right)\right)\lor \left(r^2 (2 r-1) \left(4 r^2-\text{crx}^2\right)\leq 0\land \left(\text{crx}^2+4 r^2\neq 0\lor \text{crx}^2 (5 r-2)+2 r \left(6 r^2-4 r+1\right)\neq 0\right)\right)\right)\right)\right)
\end{equation*} crx 2 + 4 ( r − 1 ) r ≤ 0 ∧ 4 r ≤ 1 ∧ ( crx = 2 r ∨ r = 0 ∨ crx + 2 r = 0 ∨ ( ( ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≤ 0 ∨ ( r 2 ( 2 r − 1 ) ( 4 r 2 − crx 2 ) ≤ 0 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ) ∨ 6 r < 1 ∨ 2 r 2 ≥ r ∨ r 2 ≤ 0 ) ∧ ( crx = 0 ∨ crx 2 r ( 4 r − 1 ) ≥ 0 ∨ crx 2 r ( 4 r − 1 ) ( crx 4 − 4 crx 2 r ( r + 1 ) + 4 r 2 ) ≥ 0 ∨ ( r 2 ( crx 10 − 4 crx 8 ( r ( 9 r − 4 ) + 1 ) + 8 crx 6 r 2 ( 62 r 2 − 48 r + 11 ) + 16 crx 4 r 2 ( 4 r ( r ( 3 r ( 3 r − 4 ) + 10 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≥ 0 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ) ∨ ( r 2 ( crx 10 + 4 crx 8 r ( 7 r − 4 ) − 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≤ 0 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 4 + 4 crx 2 r ( 3 r − 2 ) + 8 r 2 ( 8 r 2 − 4 r + 1 ) = 0 ) ∧ ( crx 10 + 4 crx 8 r ( 7 r − 4 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) + 64 r 6 = 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) ∨ crx 2 r 2 + 4 r 4 ≤ 0 ∨ ( crx 2 + 4 r 2 = 0 ∧ crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 = 0 ) ) ) ∨ ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 = 0 ∧ crx 2 r ( 4 r − 1 ) ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 ) ≤ 0 ) ∨ ( crx 2 r ( 4 r − 1 ) ( crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ) ≤ 0 ∧ ( r 2 ≤ 0 ∨ crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 = 0 ) ) ) ∧ ( crx = 0 ∨ ( crx 2 = 2 r ∧ r > 0 ) ∨ ( ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ∧ crx 8 r 2 + 16 crx 2 r 6 ≤ 4 crx 4 ( crx 2 + 1 ) r 4 ) ∨ r = 0 ∨ ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 = 0 ∧ crx 2 r 2 ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ) ≥ 0 ∧ crx 2 r ( 6 r − 1 ) ≤ 0 ) ∨ crx 2 r ( 2 r − 1 ) ≥ 0 ∨ crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ) ∧ ( 2 r 2 ≥ r ∨ ( crx 2 − 2 r ) 2 ≤ 0 ∨ r ( 2 r − 1 ) ( crx 2 − 2 r ) > 0 ) ∧ ( ( r > 0 ∧ ( r ( 2 r − 1 ) ( crx 2 − 2 r ) < 0 ∨ ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) > 0 ∧ 6 r < 1 ) ) ) ∨ crx 2 + 4 r 2 = 0 ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) ≤ 0 ∨ crx 8 r 2 + 64 r 8 ≤ 4 crx 4 ( 4 r 6 + r 4 ) ) ) ) ∧ ( crx = 2 r ∨ r = 0 ∨ crx + 2 r = 0 ∨ ( ( crx = 0 ∨ crx 2 r ( 4 r − 1 ) ≥ 0 ∨ crx 2 r ( 4 r − 1 ) ( crx 4 − 4 crx 2 r ( r + 1 ) + 4 r 2 ) ≥ 0 ∨ ( r 2 ( crx 10 − 4 crx 8 ( r ( 9 r − 4 ) + 1 ) + 8 crx 6 r 2 ( 62 r 2 − 48 r + 11 ) + 16 crx 4 r 2 ( 4 r ( r ( 3 r ( 3 r − 4 ) + 10 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≥ 0 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ) ∨ ( r 2 ( crx 10 + 4 crx 8 r ( 7 r − 4 ) − 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) − 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) + 64 r 6 ) ≤ 0 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 4 + 4 crx 2 r ( 3 r − 2 ) + 8 r 2 ( 8 r 2 − 4 r + 1 ) = 0 ) ∧ ( crx 10 + 4 crx 8 r ( 7 r − 4 ) + 16 crx 4 r 2 ( 4 r ( r ( 5 r ( 5 r − 4 ) + 11 ) − 3 ) + 1 ) + 64 r 6 = 8 crx 6 r 2 ( 2 r ( r + 8 ) − 7 ) + 16 crx 2 r 4 ( 8 r ( 5 r − 2 ) + 3 ) ∨ crx 2 r 2 + 4 r 4 ≤ 0 ∨ ( crx 2 + 4 r 2 = 0 ∧ crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 = 0 ) ) ) ∨ ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 = 0 ∧ crx 2 r ( 4 r − 1 ) ( crx 4 + 4 crx 2 r ( 7 r − 3 ) + 4 r 2 ) ≤ 0 ) ∨ crx 2 r ( 4 r − 1 ) ( crx 4 + crx 2 ( 4 ( 2 − 5 r ) r − 2 ) + 4 r 2 ) ≥ 0 ) ∧ ( crx = 0 ∨ ( crx 2 = 2 r ∧ r > 0 ) ∨ ( ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ∧ crx 8 r 2 + 16 crx 2 r 6 ≤ 4 crx 4 ( crx 2 + 1 ) r 4 ) ∨ ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 = 0 ∧ crx 2 r 2 ( crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ) ≥ 0 ∧ crx 2 r ( 6 r − 1 ) ≤ 0 ) ∨ crx 2 r ( 2 r − 1 ) ≥ 0 ∨ crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ∨ crx 2 r 2 > 0 ) ∧ ( 6 r 2 ≤ r ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≤ 0 ∨ r ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) ) > 0 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) < 0 ) ∧ ( ( r > 0 ∧ ( r ( 2 r − 1 ) ( crx 2 − 2 r ) < 0 ∨ ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) > 0 ∧ 6 r < 1 ) ) ) ∨ crx 2 + 4 r 2 = 0 ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) ≤ 0 ∨ crx 2 r 2 + 4 r 4 > 0 ∨ crx 8 r 2 + 64 r 8 ≤ 4 crx 4 ( 4 r 6 + r 4 ) ) ) ) ∧ ( r = 0 ∨ crx 2 ≥ 4 r 2 ∨ ( ( crx = 0 ∨ crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ∨ ( ( crx 2 r ( 6 r − 1 ) ≤ 0 ∨ ( r > 0 ∧ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≤ 0 ) ) ∧ ( crx 2 r ( 6 r − 1 ) < 0 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 = 0 ) ) ∨ crx 2 r ( 2 r − 1 ) ≥ 0 ∨ ( crx 8 r 2 + 16 crx 2 r 6 ≤ 4 crx 4 ( crx 2 + 1 ) r 4 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ∧ ( crx 2 ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ( crx 2 + 4 r 2 ) 2 ≥ 0 ∨ ( crx = 2 r ∧ crx 2 = 2 r ∧ crx + 2 r = 0 ∧ crx 2 + 2 r = 0 ) ) ) ) ∧ ( 6 r = 1 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 4 r 2 ≥ 0 ∨ r ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) ) > 0 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) < 0 ) ∧ ( crx = 2 r ∨ ( 0 < r < 6 1 ∧ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) > 0 ) ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ ( crx 2 − 1 ) r 2 ( crx 4 − 16 r 4 ) ≥ 0 ) ∧ ( crx = 2 r ∨ 0 < r < 6 1 ∨ crx 4 + 4 crx 2 r ( 8 r − 3 ) + 8 r 2 ( 6 r ( 3 r − 1 ) + 1 ) < 0 ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ ( crx 2 − 1 ) r 2 ( crx 4 − 16 r 4 ) ≥ 0 ) ∧ ( crx = 2 r ∨ r ( crx 2 ( 10 r − 3 ) + 2 r ( 2 r ( 8 r − 3 ) + 1 ) ) < 0 ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ 0 < r < 6 1 ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) < 0 ∨ 4 crx 4 ( 4 r 6 + r 4 ) ≤ crx 8 r 2 + 64 r 8 ) ∧ ( crx = 2 r ∨ r ( 2 r − 1 ) ( crx 2 − 2 r ) < 0 ∨ crx + 2 r = 0 ∨ crx 2 + 4 r 2 = 0 ∨ ( crx 2 ( 3 − 10 r ) < 2 r ( 2 r ( 8 r − 3 ) + 1 ) ∧ 6 r < 1 ) ∨ r 2 ( 2 r − 1 ) ( 16 r 4 − crx 4 ) < 0 ∨ 4 crx 4 ( 4 r 6 + r 4 ) ≤ crx 8 r 2 + 64 r 8 ) ∧ ( 6 r < 1 ∨ 2 r 2 ≥ r ∨ ( crx 2 − 1 ) r 2 ( crx 2 − 4 r 2 ) ≥ 0 ∨ ( ( crx 2 + 4 r 2 = 0 ∨ crx 4 + 4 crx 2 r ( 3 r − 2 ) + 8 r 2 ( 8 r 2 − 4 r + 1 ) = 0 ) ∧ ( ( crx = 2 r ∧ crx + 2 r = 0 ) ∨ crx 2 r 2 + 4 r 4 ≤ 0 ) ) ∨ ( r 2 ( 2 r − 1 ) ( 4 r 2 − crx 2 ) ≤ 0 ∧ ( crx 2 + 4 r 2 = 0 ∨ crx 2 ( 5 r − 2 ) + 2 r ( 6 r 2 − 4 r + 1 ) = 0 ) ) ) ) ) ( 22 )
若在命题 (21) 成立的条件下, 进一步让 (18) 加入大前提, 则期待 (16) 能保证 (11). 用 MMA 验证,
Resolve[ ForAll[ { cx, sx, cy} , cx^ 2 + sx^ 2 == 1 & & - 1 <= cy <= 1 , ( crx* sx) ^ 2 + ( 2 * cx* r + 2 * cy* r - 2 * r - 2 * r + 1 ) ^ 2 <= 1 ] , Reals]
FullSimplify[ %, 0 < Abs [ crx] <= 2 r & & Abs [ - 2 * r - 2 * r + 1 ] + 2 r + 2 r <= 1 ]
得
Cr x 2 + 36 r 2 ≤ 12 r . (23) \tag{23}
\text{Cr}_x^2 + 36 r^2 \leq 12 r. Cr x 2 + 36 r 2 ≤ 12 r . ( 23 )
即 MMA 认为, 在 (16)(18)(21) 成立的大前提下, 命题 (11) 等价于命题 (23), 而后者并非恒成立; 这就否定了 (16) 对 (11) 的充分性, 表明 MMA 不支持本节的主要结果 (16). 那么, 问题出在哪? 请求读者帮助! 文末附有联系方式 .
Problem 1.3
本节讨论格式 (3) 的精度. 格式 (3a) 的 modified equation
B t + u ⋅ ∇ B − ∇ ⋅ ( D ∇ B ) + K r B = − τ 2 B t 2 − h x 2 6 u B x 3 − h y 2 6 v B y 3 + h x 2 12 D x B x 4 + h y 2 12 D y B y 4 + O ( τ 2 + h x 2 + h y 2 ) . \begin{equation*}
\tag{24}
\begin{aligned}
& \qquad B_t + \boldsymbol{u} \cdot \nabla B - \nabla \cdot (D \nabla B) + K_{\text{r}} B \\
& = -\frac{\tau}{2} B_{t^2} - \frac{h_x^2}{6} u B_{x^3} - \frac{h_y^2}{6} v B_{y^3} + \frac{h_x^2}{12} D_x B_{x^4} + \frac{h_y^2}{12} D_y B_{y^4} + \mathcal{O}(\tau^2 + h_x^2 + h_y^2).
\end{aligned}
\end{equation*} B t + u ⋅ ∇ B − ∇ ⋅ ( D ∇ B ) + K r B = − 2 τ B t 2 − 6 h x 2 u B x 3 − 6 h y 2 v B y 3 + 12 h x 2 D x B x 4 + 12 h y 2 D y B y 4 + O ( τ 2 + h x 2 + h y 2 ) . ( 24 )
试图把 (24) 右端的时间误差项主部系数 B t 2 B_{t^2} B t 2 ∂ t \partial_t ∂ t ∇ \nabla ∇
B t 2 + u ⋅ ∇ B t − ∇ ⋅ ( D ∇ B t ) + K r B t = O ( τ + h x 2 + h y 2 ) , \begin{equation*}
\tag{25}
B_{t^2} + \boldsymbol{u} \cdot \nabla B_t - \nabla \cdot (D \nabla B_t) + K_{\mathrm{r}} B_t = \mathcal{O}(\tau + h_x^2 + h_y^2),
\end{equation*} B t 2 + u ⋅ ∇ B t − ∇ ⋅ ( D ∇ B t ) + K r B t = O ( τ + h x 2 + h y 2 ) , ( 25 )
∇ B t + ( u ⋅ ∇ ) ∇ B − ∇ [ ∇ ⋅ ( D ∇ B ) ] + K r ∇ B = O ( τ + h x 2 + h y 2 ) . \begin{equation*}
\tag{26}
\nabla B_{t} + (\boldsymbol{u} \cdot \nabla) \nabla B - \nabla \left[ \nabla \cdot (D \nabla B) \right] + K_{\mathrm{r}} \nabla B = \mathcal{O}(\tau + h_x^2 + h_y^2).
\end{equation*} ∇ B t + ( u ⋅ ∇ ) ∇ B − ∇ [ ∇ ⋅ ( D ∇ B ) ] + K r ∇ B = O ( τ + h x 2 + h y 2 ) . ( 26 )
由 (25)(26) 和 (24) 得
B t 2 = ( u ⋅ ∇ ) [ ( ( u ⋅ ∇ ) − K r ) B − ∇ ⋅ ( D ∇ B ) ] − K r [ u ⋅ ∇ B − ∇ ⋅ ( D ∇ B ) + K r B ] − ∇ ⋅ [ D ( u ⋅ ∇ ) ∇ B − D ∇ [ ∇ ⋅ ( D ∇ B ) ] + K r D ∇ B ] + O ( τ + h x 2 + h y 2 ) . \begin{equation*}
\tag{27}
\begin{aligned}
B_{t^2} & = (\boldsymbol{u} \cdot \nabla) \left[ \left( (\boldsymbol{u} \cdot \nabla) - K_{\text{r}} \right) B - \nabla \cdot (D \nabla B) \right] \\
& - K_{\mathrm{r}} \left[ \boldsymbol{u} \cdot \nabla B - \nabla \cdot (D \nabla B) + K_{\text{r}} B \right] \\
& - \nabla \cdot \left[ D (\boldsymbol{u} \cdot \nabla) \nabla B - D \nabla \left[ \nabla \cdot (D \nabla B) \right] + K_{\mathrm{r}} D \nabla B \right] \\
& + \mathcal{O}(\tau + h_x^2 + h_y^2).
\end{aligned}
\end{equation*} B t 2 = ( u ⋅ ∇ ) [ ( ( u ⋅ ∇ ) − K r ) B − ∇ ⋅ ( D ∇ B ) ] − K r [ u ⋅ ∇ B − ∇ ⋅ ( D ∇ B ) + K r B ] − ∇ ⋅ [ D ( u ⋅ ∇ ) ∇ B − D ∇ [ ∇ ⋅ ( D ∇ B ) ] + K r D ∇ B ] + O ( τ + h x 2 + h y 2 ) . ( 27 )
Eq. (27) 有些复杂, 感兴趣的读者可尝试将其代入 (24), 整理出只含空间导数的截断误差项主部, 看有无可能通过适当选取网格参数 τ , h x , h y \tau, h_x, h_y τ , h x , h y O ( τ ) \mathcal{O}(\tau) O ( τ ) O ( τ + h x 2 + h y 2 ) \mathcal{O}(\tau + h_x^2 + h_y^2) O ( τ + h x 2 + h y 2 )
Problem 1.4
将格式 (3) 写成迭代形式,
[ B O ] i , j n + 1 = ( 1 − 2 r x − 2 r y ) [ B O ] i , j n − τ [ K r K r K a ] [ B O ] i , j n + τ [ q i , j K a O s a t ] + ∑ k ∈ { − 1 , 1 } [ ( r x − C r x 2 s g n k ) [ B O ] i + k , j n + ( r y − C r y 2 s g n k ) [ B O ] i , j + k n ] , \begin{equation*}
\tag{28}
\begin{aligned}
\begin{bmatrix*}
B \\
O
\end{bmatrix*}^{n+1}_{i,j} & = (1 - 2 r_x - 2 r_y) \begin{bmatrix*}
B\\
O
\end{bmatrix*}^n_{i,j} - \tau \begin{bmatrix*}
K_{\mathrm{r}} & \\
K_{\mathrm{r}} & K_{\mathrm{a}}
\end{bmatrix*} \begin{bmatrix*}
B \\
O
\end{bmatrix*}^{n}_{i,j} + \tau \begin{bmatrix}
q_{i,j} \\
K_{\mathrm{a}} O_{\mathrm{sat}}
\end{bmatrix} \\
& + \sum_{k \in \{ -1, 1 \}} \left[
\left( r_x - \frac{\mathrm{Cr}_x}{2} \mathop{\mathrm{sgn}}{k} \right) \begin{bmatrix}
B \\
O
\end{bmatrix}^n_{i+k, j} + \left( r_y - \frac{\mathrm{Cr}_y}{2} \mathop{\mathrm{sgn}}{k} \right) \begin{bmatrix}
B \\
O
\end{bmatrix}^n_{i, j+k}
\right],
\end{aligned}
\end{equation*} [ B O ] i , j n + 1 = ( 1 − 2 r x − 2 r y ) [ B O ] i , j n − τ [ K r K r K a ] [ B O ] i , j n + τ [ q i , j K a O sat ] + k ∈ { − 1 , 1 } ∑ [ ( r x − 2 Cr x sgn k ) [ B O ] i + k , j n + ( r y − 2 Cr y sgn k ) [ B O ] i , j + k n ] , ( 28 )
其中 BOD 的源项
q i , j = ∑ ( x 0 , y 0 ) Q B ( x 0 , y 0 ) h x h y I i , j ( x 0 , y 0 ) , \begin{equation*}
\tag{29}
q_{i,j} = \sum_{(x_0, y_0)} \frac{Q_{\mathrm{B}}(x_0, y_0)}{h_x h_y} \mathrm{I}_{i,j} (x_0, y_0),
\end{equation*} q i , j = ( x 0 , y 0 ) ∑ h x h y Q B ( x 0 , y 0 ) I i , j ( x 0 , y 0 ) , ( 29 )
where ( x 0 , y 0 ) (x_0, y_0) ( x 0 , y 0 ) Q B ( x 0 , y 0 ) Q_{\mathrm{B}}(x_0, y_0) Q B ( x 0 , y 0 ) I i , j ( x 0 , y 0 ) \mathrm{I}_{i, j}(x_0, y_0) I i , j ( x 0 , y 0 ) Matlab 实现数值求解二维 Streeter-Phelps 模型 (1)(2) 的 FTCS 方案 (28).
在 (28) 中, 令
[ B O ] n + 1 = [ B O ] n = : [ B ∗ O ∗ ] , \begin{equation*}
\tag{30}
\begin{bmatrix*}
B \\
O
\end{bmatrix*}^{n+1} = \begin{bmatrix*}
B \\
O
\end{bmatrix*}^{n} =: \begin{bmatrix*}
B^{*} \\
O^{*}
\end{bmatrix*},
\end{equation*} [ B O ] n + 1 = [ B O ] n =: [ B ∗ O ∗ ] , ( 30 )
得数值稳态 (steady-state) 解 B ∗ , O ∗ B^*, O^* B ∗ , O ∗
τ [ q i , j K a O s a t ] = ( 2 r x + 2 r y ) [ B ∗ O ∗ ] i , j + τ [ K r K r K a ] [ B ∗ O ∗ ] i , j − ∑ k ∈ { − 1 , 1 } [ ( r x − C r x 2 s g n k ) [ B ∗ O ∗ ] i + k , j + ( r y − C r y 2 s g n k ) [ B ∗ O ∗ ] i , j + k ] . \begin{equation*}
\tag{31}
\begin{aligned}
\tau \begin{bmatrix}
q_{i,j} \\
K_{\mathrm{a}} O_{\mathrm{sat}}
\end{bmatrix} & = (2 r_x + 2 r_y) \begin{bmatrix*}
B^{*} \\
O^{*}
\end{bmatrix*}_{i,j} + \tau \begin{bmatrix*}
K_{\mathrm{r}} & \\
K_{\mathrm{r}} & K_{\mathrm{a}}
\end{bmatrix*} \begin{bmatrix*}
B^{*} \\
O^{*}
\end{bmatrix*}_{i,j} \\
& - \sum_{k \in \{ -1, 1 \}} \left[
\left( r_x - \frac{\mathrm{Cr}_x}{2} \mathop{\mathrm{sgn}}{k} \right) \begin{bmatrix}
B^{*} \\
O^{*}
\end{bmatrix}_{i+k, j} + \left( r_y - \frac{\mathrm{Cr}_y}{2} \mathop{\mathrm{sgn}}{k} \right) \begin{bmatrix}
B^{*} \\
O^{*}
\end{bmatrix}_{i, j+k}
\right].
\end{aligned}
\end{equation*} τ [ q i , j K a O sat ] = ( 2 r x + 2 r y ) [ B ∗ O ∗ ] i , j + τ [ K r K r K a ] [ B ∗ O ∗ ] i , j − k ∈ { − 1 , 1 } ∑ [ ( r x − 2 Cr x sgn k ) [ B ∗ O ∗ ] i + k , j + ( r y − 2 Cr y sgn k ) [ B ∗ O ∗ ] i , j + k ] . ( 31 )
要求解 (31), 可先求解 (31) 中关于 { B i , j ∗ } \{ B^*_{i,j} \} { B i , j ∗ } sparse()B ∗ B^* B ∗ { O i , j ∗ } \{ O^*_{i,j} \} { O i , j ∗ }
[ B O ] n → [ B ∗ O ∗ ] , ( n → ∞ ) . \begin{equation*}
\tag{32}
\begin{bmatrix*}
B \\
O
\end{bmatrix*}^n \to \begin{bmatrix*}
B^{*} \\
O^{*}
\end{bmatrix*}, \quad (n \to \infty).
\end{equation*} [ B O ] n → [ B ∗ O ∗ ] , ( n → ∞ ) . ( 32 )
为确保 FTCS 格式 (28) 稳定且收敛, 应适当选取网格参数 h x , h y h_x, h_y h x , h y τ \tau τ 1.3 节的分析, 当 v = 0 , D x = D y = D v = 0, \, D_x = D_y = D v = 0 , D x = D y = D
h x = h y = : h ≤ 2 D / ∣ u ∣ , τ ≤ ( 4 D / h 2 + max { K a , K r } ) − 1 . \begin{equation*}
\tag{33}
h_x = h_y =: h \le {2 D / |u|}, \quad \tau \le \left(4 D / h^2 + \max{\{ K_{\mathrm{a}}, K_{\mathrm{r}} \}} \right)^{-1}.
\end{equation*} h x = h y =: h ≤ 2 D /∣ u ∣ , τ ≤ ( 4 D / h 2 + max { K a , K r } ) − 1 . ( 33 )
由 (33), 我们选择 h = 5 × 1 0 − 2 m , τ = 1.24 × 1 0 − 3 s h = 5 \times 10^{-2} \, \text{m}, \, \tau = 1.24 \times 10^{-3} \, \text{s} h = 5 × 1 0 − 2 m , τ = 1.24 × 1 0 − 3 s
我们测试了以下四种排放方案:
单点排放 (single-point discharge). 所有 BOD 都排放在 ( x , y ) = ( 10 , 15 ) m (x, y) = (10, 15) \, \text{m} ( x , y ) = ( 10 , 15 ) m 两点排放 (two-points discharge). BOD 被平均地排放在 ( x , y ) = ( 10 , 15 ) , ( 70 , 15 ) (x, y) = (10, 15), (70, 15) ( x , y ) = ( 10 , 15 ) , ( 70 , 15 ) 线分布排放 (line-distributed discharge). BOD 排放源均匀分布在线段 x ∈ [ 10 , 70 ] , y = 15 x \in [10, 70], \, y = 15 x ∈ [ 10 , 70 ] , y = 15 面分布排放 (face-distributed discharge). BOD 排放源均匀分布在区域 ( x , y ) ∈ [ 10 , 70 ] × [ 10 , 20 ] (x, y) \in [10, 70] \times [10, 20] ( x , y ) ∈ [ 10 , 70 ] × [ 10 , 20 ]
测试结果示于下面 8 个图. 每个图都标注了相应时刻的 BOD 或 DO 的最值. 发现以下特点:
对于点排放 (单点, 两点) 情形, BOD 最高值逐渐升高至稳态值, 且位于最下游的排放源点. DO 最低值逐渐下降至稳态值, 且位置随时间推移有向下游移动的趋势, 最终稳定在排放源点的下游.
对于均匀分布排放 (线, 面) 情形, BOD 最高值逐渐升高至稳态值, 且位置随时间推移有从源区中心向下游移动的趋势, 最终稳定在源区下游边界的内侧. DO 最低值逐渐下降至稳态值. 随时间推移, DO 最低值点有从源区中心向下游移动的趋势, 且相对于 BOD 最高值点, 有从滞后 (即位于相对上游) 到超前的移动趋势.
采取分布排放, 可以大幅降低 BOD 峰值, 因为相比点排放可以大幅降低排放源强度. 但未必能使 DO 最低值达到最高. 物理解释: 点排放引起 DO 在小范围内有快速下降的倾向, 分布排放引起 DO 在大范围内有缓慢下降的倾向. 若平流速度充分大或点排放源强度相对足够小, 则对于点排放情形, 水团受点排放源的影响时间较小, 引起的 DO 损失足够被水-气氧气通量来补偿; 而对于分布排放, 水团在较长时间内持续地被 BOD 反应掉, 足以产生比点排放情形更低的 DO 最低值.
